package org.cosmocalc;

/**
 * @author michaelw
 * 
 */
public class CosmoCalculator {
	public Cosmology cosmology;
	public Drdz drdzObj;
	public Drdzprop drdzpropObj;
	public double ageOfU; 
	
	public CosmoCalculator() {
		cosmology = new Cosmology();
		drdzObj = new Drdz(cosmology);
		drdzpropObj = new Drdzprop(cosmology);
		ageOfU = this.time(0.);
	}

	public void setCosmo(double om, double ol, double hub) {
		cosmology.setCosmo(om, ol, hub);
		ageOfU = this.time(0.);
	}
	/**
	 * Calculate the age of the Universe at a given redshift.
	 * @param z redshift
	 * @return Age in Giga-years.
	 */
	public double time(double z) {
		double c = 2.99792458E10 / 3.08568E24 * 3.155815E7;
		// (cm/s)/(3.08568e24cm/Mpc)*(3.16e7s/yr) = Mpc/yr
		if (z < 250.) {
			return (-1. / c * Simpson.qsimp(drdzpropObj, 250., z) / 1.e9); // Giga-years
		}
		return (0.);
	}
	/**
	 * Calculate the lookback time at a given redshift.
	 * @param z redshift
	 * @return Lookback time in Giga-years.
	 */
	public double lookbackTime(double z) {
		return (ageOfU - time(z));
	}

	/**
	 * Calculate the angular distance at a given redshift.
	 * @param z redshift
	 * @return Angular distance in kpc per arcsec.
	 */
	public double angDist(double z) {
		return ((D(z) / (1. + z))/206265.*1000.);
	}
	/**
	 * Calculate the luminosity distance at a given redshift.
	 * @param z redshift
	 * @return Luminosity distance in Gpc.
	 */
	public double lumDist(double z) {
		return (D(z) * (1. + z))/1.E3;
	}

	/**
	 * Calculate volume at a given redshift.
	 * @param z redshift
	 * @return
	 */
	public double dV(double z) {
		return (Math.pow(D(z), 2.) * drdzObj.f(z));
	}

	private double r(double z) {
		return (Simpson.qsimp(drdzObj, 0., z));
	}

	private double D(double z) {
		double Rtilt;
		if (cosmology.getKappa() < 0.) {
			Rtilt = 1. / Math.sqrt(-cosmology.getKappa());
			return (Rtilt * sinh(r(z) / Rtilt));
		}
		if (cosmology.getKappa() > 0.) {
			Rtilt = 1. / Math.sqrt(cosmology.getKappa());
			System.out.println(Rtilt);
			return (Rtilt * Math.sin(r(z) / Rtilt));
		}
		if (cosmology.getKappa() == 0.) {
			return (r(z));
		}
		return -1.;
	}

	private double sinh(double x) {
		return ((Math.exp(x) - Math.exp(-x)) / 2.);
	}

}
